#ifndef _PACKAGE_
#define _PACKAGE_
#include "../DF_Common/basis_state.h"
#include "../DF_Common/com_sparse_row.h"
#include "../DF_Common/coupling_coeff.h"
#include "../DF_Common/progress_bar.h"
#include "../DF_Common/const_nucl.h"
#include "../DF_Local/antisym_coeff.h"
#include "../DF_Local/cont_term.h"
#include "../DF_Local/ope_term.h"
#include "../DF_Local/tpe_term.h"
#include "../DF_Local/trans_coeff.h"

#include <cmath>
#include <fstream>
#include <iostream>
#include <string>
using std::cin;
using std::cout;
using std::endl;
using std::string;
using std::vector;

class Package {
public:
  Package(TripleP &TP_x, Jacobi_NL &jacNL_x, Jacobi &jac_x, TransCoeff &trans_x,
          Antisym &anti_x, Cont_term &cont_x, Ope_term &ope_x, Tpe_term &tpe_x)
      : TP(TP_x), jacNL(jacNL_x), jac(jac_x), trans(trans_x), anti(anti_x),
        cont(cont_x), ope(ope_x), tpe(tpe_x) {}
  Package(int hw_t, TripleP &TP_x, Jacobi_NL &jacNL_x, Jacobi &jac_x,
          TransCoeff &trans_x, Antisym &anti_x, Cont_term &cont_x,
          Ope_term &ope_x, Tpe_term &tpe_x)
      : Hbar_Omega(hw_t), TP(TP_x), jacNL(jacNL_x), jac(jac_x), trans(trans_x),
        anti(anti_x), cont(cont_x), ope(ope_x), tpe(tpe_x) {}
  TripleP &TP;
  Jacobi_NL &jacNL;
  Jacobi &jac;
  TransCoeff &trans;
  Antisym &anti;
  Cont_term &cont;
  Ope_term &ope;
  Tpe_term &tpe;
  int Hbar_Omega;

  double LEC_C1_P;
  double LEC_C3_P;
  double LEC_C4_P;
  double LEC_CD_P;
  double LEC_CE_P;
  double LEC_gA_P;

  SMatrix_Ch vec_S, vec_Tz_S;
  SMatrix_Ch vec_flag_S, vec_flag_Tz_S;       // sparse matrix for vec
  vector<vector<vector<double>>> vec;         // vec[TP_channel][i][j]
  vector<vector<vector<double>>> vec_Tz;      // vec[TP_channel_Tz][i][j]
  vector<vector<double>> vec_TzM;             // vec[i][j]
  vector<vector<vector<double>>> JacA_vec;    // JacA_vec[jaco_channel][i][j]
  vector<vector<vector<double>>> vec_anti;    // vec[TP_channel][i][j]
  vector<vector<vector<double>>> vec_anti_Tz; // vec[TP_channel][i][j]

  vector<vector<vector<int>>> vec_flag;    // 0: not calculate 1: calculated
  vector<vector<vector<int>>> vec_flag_Tz; // 0: not calculate 1: calculated

  double cal(State_TP &TP_f, State_TP &TP_i);
  double cal2(State_TP &TP_f, State_TP &TP_i);

  void build(int TP_ch);
  void build();
  void build_JacAvec();
  void print_JacAvec(std::string filename);
  void print(std::string filename);
  void print_Tz(std::string filename);
  void print_TzM(std::string filename);

  double Check_Anti_T(SixJ_coeff &sixJ, SingleP &SP, int J2, int T2, int flag,
                      int index_f, int index_i);
  double Check_Anti_T_2(SixJ_coeff &sixJ, SingleP &SP, int J2, int T2, int flag,
                        int index_f, int index_i);
  double Check_Anti_Tz(SixJ_coeff &sixJ, SingleP &SP, int J2, int T_2z,
                       int flag, int index_f, int index_i);
  // flag = 0 : a<->b; flag = 1: b <-> c
  void build_Tz(CG_coeff &cg, SingleP &SP);
  void build_S_Tz(CG_coeff &cg, SingleP &SP, SMatrix_Ch &vec_in,
                  SMatrix_Ch &vec_flag_in, SMatrix_Ch &vec_out,
                  SMatrix_Ch &vec_flag_out);
  void build_TzM(CG_coeff &cg, SingleP &SP);
  int FindTPx(int ch, int T2ab, int J2ab, int T2z, int a, int b, int c);
  int FindTPx_Tz(int J2, int T_2z, int J2ab, int a, int b, int c);

  void build_anti_bc(SixJ_coeff &sixJ, SingleP &SP);
  void build_anti_ab_Tz(SixJ_coeff &sixJ, SingleP &SP);
  void build_anti_bc_Tz(SixJ_coeff &sixJ, SingleP &SP);

  //
  vector<State_SP> sp_nnlosat_state;
  vector<int> z2sat;
  vector<int> sat2z;
  vector<vector<vector<double>>> nnlosat_vec; // vec[TP_channel][i][j]
  void print_nnlosat(SingleP SP_t, std::string sp_nnlosat_file,
                     std::string filename);
  void Readnnlosat(SingleP SP_t, std::string sp_nnlosat_file,
                   std::string me_nnlosat_file, vector<int> abc_range);
  void Readnnlosat_SP(SingleP SP_t, std::string sp_nnlosat_file);
  int FindTP_index(int a, int b, int c, int J2ab, int T2ab, int J2, int T2);
  void CompareVsNNlosat(string filename, string filename2, string filename3);
  void Compare_satVSzhuo(string filename, string filename2,
                         vector<int> abc_rang);

private:
  int phase_cal(int n) { return n % 2 ? -1 : 1; }
  //
};

#endif
